Impact of high-energy tails on granular gas properties 
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The velocity distribution function of granular gases in the homogeneous cooling state as well as 
some heated granular gases decays for large velocities as / oc exp(— const. u). That is, its high-energy 
tail is overpopulated as compared with the Maxwell distribution. At the present time, there is no 
theory to describe the influence of the tail on the kinetic characteristics of granular gases. We develop 
an approach to quantify the overpopulated tail and analyze its impact on granular gas properties, 
in particular on the cooling coefficient. We observe and explain anomalously slow relaxation of the 
velocity distribution function to its steady state. 

PACS numbers: 45.70.-n,51.10.+y 



When a homogeneous granular gas evolves in the ab- 
sence of external forces, it develops a velocity distribution 
similar to a molecular gas, however, its temperature de- 
cays due to the dissipative nature of particle collisions. 
The velocity distribution function of granular gases has 
attracted much scientific attention since it deviates char- 
acteristically from the Maxwell distribution. There are 
two types of deviations: First, there are deviations from 
the Maxwellian in the main part of the distribution [l|, Q j 
where the particle velocities are close to the thermal 
velocity vrit) — \/2T (unit particle mass is assumed). 
Second, the high-energy part, v ^ vt, deviates in its 
functional form, i.e., the distribution function decays as 
/ oc exp(— const. u) 0, instead of / cx exp(— const. u^) 
as expected for the Maxwell distribution. Heated gran- 
ular gases with a Gaussian thermostat are equivalent to 
gases in the homogeneous cooling state (HCS) That 
is, the addressed properties apply for a wide class of gran- 
ular gases. Figure [T] shows the velocity distribution func- 
tion and its deviation from the Maxwellian (details see 
below). 

Both types of deviations are characterized by the co- 
efficient of restitution e describing the post-collision par- 
ticle velocities vl and V2 as functions of the pre-coUision 
velocities, 



H/2 = ^1/2 T 
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— V2 and the unit vector e 
at the moment of the collision. 



with V12 = Vl 
in ~ r2) I |ri - r2\ 

In the HCS, the velocity distribution f {v,t) (where 
the time r is measured in the average number of colli- 
sions per particle) can be reduced to a time-independent 
distribution function / (c ) by the transformation 
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FIG. 1: Velocity distribution function of a granular gas / (c) 
as normal and logarithmic plot. The symbols show a simu- 
lation of A'' = 10* particles for e = 0.3. For c ~ 1, / (c) is 
well described by the second order Sonine expansion, Eq. Q 
(top). For comparison the Maxwell distribution is also shown. 
For c ^ 1 the distribution function decays exponentially slow 
(bottom), see Eq. (|6]). The tail starts at c « c*, Eq. (|9]). 



with Haff 's law [5i] for the temperature evolution, 
dT 

— = -27T , i.e. T(t) = r(0) exp(-27r) . (3) 
dr 

The main part of the distribution function, c ~ 1, can 
be described with good accuracy by a second order So- 
nine polynomials expansion around the Maxwell distri- 
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bution (j){c) = TT ■^/^ exp(--c^): 

/(c) = <^(c) [1 + ai5i (c2) + 02^2 (c") + 



(4) 



It can be shown that ai — 0, therefore, the leading devia- 
tions from the Maxwell distribution are due to the second 
Sonine polynomial <5'2(c^) = cVS — 5c^/2 + 15/8 and the 
respective coefficient a2 P, 0, @| : 



a2 = 



16(l-£)(l-2£2) 
81- 17£ + 30e2(i-e) 



(5) 



The good agreement of Eq. (|4]) with simulation data in 
the region c ~ 1 can be seen in Fig. [1] (top) . 

For c ^ 1, Eq. ^ fails to represent the velocity dis- 
tribution function due to its different functional form. It 
has been shown that for particles with high velocities the 
distribution function develops an exponential tail 0,0] j 
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which is illustrated in the bottom part of Fig. (TJ here /i2 
is the second moment of the collision integral, e.g. p]- 

The overpopulation of the tail is a rather general fea- 
ture of granular gases: After theoretically predicted^ 
it was found for gases in the HCS also numerically 
as well as for driv en g ases and was also detected experi- 
mentally, e.g. 1^, 11, H, 13, 14 1 . In spite of its obvious 



importance, still it lacks a theoretical description which 
allows to quantify its impact on granular gas properties. 
Indeed, neither the numerical prefactor B in Eq. ^ is 
known, nor the threshold velocity above which the tail 
is overpopulated. In the present letter we address this 
problem numerically and analytically. We develop an 
approach to quantify the high-energy tail and estimate 
its impact on gas properties. The impact of the tail on 
the cooling rate is studied in detail. 

\Ve perform Direct Simulation Monte Carlo (DSMC) 
H [la m of iV = 10^ granular particles. DSMC is par- 
ticularly suited to simulate large systems over long time 
in the HCS, that is, to suppress spatial correlations which 
give rise to vortices [3| and clusters [l3|- We started 
at r(0) = 1 and simulated until the particle velocities 
approached the double precision number representation, 
i.e., until T k, 10"^^. For e = 0.9 this corresponds to a 
total of 5 X 10^° collisions or 1,000 collisions per parti- 
cle. Neglecting the first 2 • 10^ collisions, after each 10® 
collisions (2 collisions per particle) we recorded a snap- 
shot of the scaled velocities, c. The distribution function, 
/(c), was then obtained by binning of up to 100 of such 
snapshots in 100 intervals. 

The numerical values of the constants B and b were 
then determined by performing a least-mean-square fit 
of the linear function B' — he, where B' = log B. To find 
the numerical value of the threshold velocity where the 
tail starts, c*, (see the discussion below) we discriminate 
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FIG. 2: The slope h of the exponential tail. For e < 0.7 we 
find good agreement of the DSMC results with the theoretical 
expression given in Eq. ([6} (full line) . For smaller dissipation 
the data deviate (see text). 



between two cases. For the first case, where the functions 
-Bexp(— 6c) and /(c) according to the Sonine expansion, 
Eqs. (|3]), (O intersect, we determine c* as the mean 
of the two intersection points; these were very close to 
each other in all our simulations. For the second case, 
where the functions do not intersect, c* was determined 
as the scaled velocity that minimizes the distance be- 
tween both functions. Such simulations were performed 
for £ = 0.1, 0.2,..., 0.9. 

In Fig. [2] the numerical results for the slope h are com- 
pared with theoretical prediction, Eq. ([6|). For £ < 0.7 
we find good agreement, whereas for larger e the data 
deviate. This can be understood from the theoretical 
argument used to derive /(c) (x exp(— 5 c) d, IIO]: For 
c ^ oo the gain term of the Boltzmann equation may 
be neglected as compared with the loss term. The larger 
the restitution coefficient £, the later the tail starts. Con- 
trary, for molecular gases where £ = 1, the gain and loss 
terms balance each other for all velocities leading to the 
Maxwell distribution, /(c) cx exp(— c^). The deviation 
for £ > 0.7 originates from the fact that the gain term 
cannot be neglected for large e for the accessible interval 
of velocities. 

To obtain the the threshold velocity c* analytically we 
assume that /(c) may be sufficiently well described by a 
combination of Eq. ([5]), valid for c ~ 1 and Eq. Q for 
the tail, 

/(c) = Ac^e-"^ [1 -f a-2.S-2.{c^)\ 9{c,-c)+Bc^e-'''e{c-c,) , 

(7) 

with the Heaviside function Q{x), i.e., we disregard the 
transient region [20| between the near-Maxwellian and 
exponential parts of the distribution. Using DSMC we 
checked the Ansatz ([7]), whose eligibility is illustrated in 
Fig. [l]for £ = 0.3. The unknown parameters A, B and 
the threshold velocity c* may be found from the nor- 
malization condition and continuity of the distribution 
function itself and its first derivative 

/(c,+ 0) = /(c,-0), /'(c,+ 0) = /'(c,-0), (8) 
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where /' = df /dc. From normalization follows then 
b a2 {2cl - 5c,) 



A- 



2 2(l + a252(c2)) 



(9) 



^3 (2 + 6c42 + 6c,))e-''^* +^Erf(c,) 
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Ak{c^) 



c,(4 + a2C^(2c2-5))e~^- (10) 

(11) 



with 



(12) 



Solving the fifth order equation ([9]) numerically for c,, 
we obtain A, k and finally B, Fig. [31 
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FIG. 3: The threshold velocity c» and the parameter B as 
obtained from DSMC (symbols) together with the solution of 
Eqs. ((9)1 and Again, for e < 0.7 we find good agreement. 

The very good agreement between simulations and the 
theoretical predictions for the coefficients A(e), B{e) and 
the transition velocity c*(e) manifests the adequacy of 
the Ansatz ([7]). The deviations for e ^ 0.7 occur already 
for the slope b of the exponential tail, Fig. [21 which is 
not related to the Ansatz. It may be explained similarly 
as the deviation of b from its theoretical value: For large 
£ > 0.7 the system of = 10^ particles is not large 
enough to develop a well-detectable exponential tail. 

Knowing the distribution function, Eq. ([7]) and its pa- 
rameters A, B and c, as functions of e, we can quantify 
the impact of the exponential tail on kinetic quantities. 
In this letter we focus on the temperature decay rate 7; 
for the diffusion coefficient, other transport coefficients 
as well as for technical details we refer to ^2^. 

The standard analysis (e.g. [l, [3|) yields for the tem- 
perature decay rate, when the stationary velocity distri- 
bution is achieved, 

7 = (i-s')T#r> (13) 
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where 



Jk^ldcidc2 /dee(-ci2-e)|ci2-e|'''+V(ci)/(c2). (14) 



Disregarding the exponential tail, the energy decay rate 
reduces to |2| 



70(e) 



1 -£21 + (3/16) 02(2) 



(15) 



6 1- (1/16) 02(2) 
We applied DSMC of 10* particles for different e and 
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FIG. 4: Inset: The temperature relaxation coefficient 7(e) 
over e as obtained by DSMC, 7dsmc, (points) and due to Eq. 
(|15p . 70. The logarithmic plot shows (70 — 7dsmc), that is, 
the inffuence of the exponential tail on the cooling coefficient. 
The dashed line shows the best exponential fit. 

recorded the temperature Tdsmc(t). Then 7dsmc was 
determined by fitting TbsMc(''') for t 1 to its asymp- 
totic law, Tdsmc oc exp(-2 7DSMC t), Eq. Figure 
m (inset) shows 7dsmc(£) (points) together with the an- 
alytical result, 7o(e), Eq. ([TCjl (fine). In this represen- 
tation we hardly see any discrepancy between theoret- 
ical and numerical data. The difference between these 
curves (main part of Fig. [3]) which quantifies the im- 
pact of the overpopulated tail on the cooling rate re- 
veals, however, a clear dependence on e. The scaling law, 
70 — 7DSMC oc exp(— 6e), shown here only as a numeri- 
cal result is, however, difficult to confirm analytically: In 
spite of the simple functional form of /(c), Eq. ([7]), an 
accurate analytical expression for the cooling coefficient 
may be obtained only in the limit e ^ 1 [1^ . As it follows 
from the discussion below, the studied system of A^ = 10* 
particles is not sufficiently large to analyze this limit by 
numerical simulations. 

So far, the discussion refers to the state when the veloc- 
ity distribution /(c) has relaxed to its stationary form. 
Now we ask the question, how fast does this happen? 
To address this problem we initialize the particle veloc- 
ities according to a Maxwell distribution at r(0) = 1 
and investigate the decay of temperature as a function of 
the average number of collisions r. Asymptotically, i.e., 
when the gas has adopted its asymptotic distribution, 
the temperature evolves according to Haff's law, Eq. ([3]). 
Thus, the time lag which is needed for a system to reach 
Haff's evolution corresponds to the relaxation time of the 
distribution function to achieve its stationary form /(c). 
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Using the coefficient 7dsmc described above, we de- 
fine the temperature 7fit(T) cx exp(— 27dsmc ■'')• By def- 
inition, for T > 1 we have Tdsmc ~ Tfit since 7dsmc 
was determined as the best exponential fit to Tdsmc (''') 
for T 1. Therefore, the quantity 1 — Tdsmc (''')/7fit(T) 
characterizes the relaxation of the distribution function 
to its stationary form. Figure [5] shows the relaxation for 
different values of the coefficient of restitution. We note 
that depending on e, the relaxation to the level of "natu- 
ral" fluctuations in the system takes approximately 20 to 
30 collisions per particle. This slow relaxation is very dif- 
ferent from that in a molecular gas, where it takes very 
few (3-5) collision per particle to develop the Maxwell 
distribution. Similar slow relaxation of the high-energy 



Using the Ansatz /(c, t) = B exp [—w{t)c] , we recast Eq. 
(fT6l) into 




FIG. 5: Relaxation of the temperature decay to Half's law for 
= 10* and e = 0.3, 0.4, 0.6, characterizing the relaxation of 
the distribution function to its stationary form. 

tail of the velocity distribution was reported for a gas of 
elastic hard spheres 21|. The relaxation mechanism for 
the case of granular gases is, however, completely differ- 
ent and depends on e: While the tail in clastic gases is 
fed by the gain term of the Boltzmann equation, for the 
dissipative gases this term is negligible. The formation 
of the tail for the latter case occurs exclusively due to 
permanent cooling of the gas, so that the scaled velocity 
of a particle, c = v/vt, increases due to decaying vt- 
By this mechanism particles enter the tail and keep it 
overpopulated. 

The relaxation time, i.e., the time of formation of the 
exponential tail, may be quantified. The relaxation to 
the stationary distribution / (c) is described by 0, [l^l 



y (3 + c^^ /(c, r) + Jo|:/(c, r) - / (/, /) (16) 

where I{f,f) is the reduced collision integral (e.g. Q) 
and Jo is defined in Eg. ((T4| . Neglecting the incoming 
term for c ^ 1 @j B Q) the collision integral may be 
approximated by 



dr 3 Jo Jo 



c > 1 



with the solution 



w{t) =b + {1 -b) exp [-T/To(e) 



(18) 



(19) 



where b — 37r//i2 coincides with Eq. ^ and T^^{e) = 
fJ-2/SJo. Neglecting 02, which characterizes small defor- 
mations of the main part of the distribution with respect 
to the Maxwellian, and the contribution from the tail, we 
obtain Jo = 2'\/27r and hence 
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(20) 



-7rc/(c) 



O 1 . 



(17) 



(see [23| for details). For e = 0.4 we obtain the relaxation 
time. To « 7.1. 

Let us compare the theoretical prediction, Eq. ((20|) . 
with the numerical results: The relaxation of (1 — 
TosMc/Tfit) as plotted in Fig. [5]reveals two stages. First, 
for r ^ 10, we observe relaxation of the main part of the 
velocity distribution, where c « 1. The initial Maxwell 
distribution relaxes here to the distribution given by the 
Sonine expansion, Eq. During the second stage the 
overpopulated tail is formed; the plotted quantity de- 
creases for e — 0.4 by a factor of 10 in the time span 
Ar — 25 ranging from r = 10 to t = 35. This leads to 
a numerical relaxation time To = 25/log(10) « 10.8, in 
agreement with the above theoretical estimates. 

The theory also predicts that the relaxation time in- 
creases with increasing e. While this tendency is con- 
firmed for e = 0.3 and e = 0.4, it is seemingly violated 
for e ~ 0.6, see Fig. [5l We argue however that this is, 
presumably, a finite size effect, which may be understood 
as follows. According to the mechanism of the tail for- 
mation, discussed above, the gain term of the collision 
integral does not contribute to the tail. Instead, parti- 
cles enter the tail due to increase of the scaled velocity, 
c= w/wt, when the thermal velocity vt decreases along 
with temperature T. The temperature decay and hence 
the formation of the tail is slower for larger e, that is, the 
relaxation time To is larger, Eq. 1201 

On the other hand, the total number of particles in the 
tail moving at velocities c > c*, decreases with increas- 
ing £. Correspondingly, the deviation of the distribution 
function from its steady state /(c), quantified here by 
(1 — Tdsmc /Tfit), becomes smaller for smaller dissipa- 
tion. Consequently, the relaxation of this quantity may 
be traced only as long as it exceeds the level of natu- 
ral fluctuations. The smaller the system, the larger is 
the impact of the fluctuations. Therefore, if the number 
of particles in the tail is not sufficiently large, the value 
of (1 — TDSMc/Tfit) drops quickly below the fluctuation 
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level, making an accurate numerical estimate of the re- 
laxation time impossible. This is the case for e = 0.6 
in Fig. [S] where a seemingly fast relaxation is observed 
due to large fluctuations. Hence, we conclude that the 
observed relaxation curves do not contradict the predic- 
tions of the theory. They indicate, however, that the 
size of the system of N = 10* particles is not sufficient 
to study the relaxation of the distribution function for 
e — 0.6 or larger. For a very large system we expect 
increase of the relaxation time tq with increasing coef- 
ficient of restitution e in agreement with the theoretical 
analysis. 

In summary, we investigated the velocity distribution 
function of a granular gas and the impact of its overpop- 
ulated high-energy tail on the cooling coefficient, which 
is the main characteristics of a granular gas in the HCS. 
We proposed a unified functional form of the distribution 
function which comprises its main part (v/vt = c ~ 1) 
whose deviation from the Maxwell distribution is de- 
scribed by the second-order Sonine expansion, and the 
overpopulated tail which decays exponentially. We de- 
rived e-dependent coefficients of the proposed Ansatz 
along with the scaled velocity c*, which separates the 
main part of the velocity distribution (c w 1) and the 
tail part (c ^ 1). For e < 0.7 the analytical results agree 
well with large scale DSMC of 10* particles, while the 
deviations for e > 0.7 may be attributed to finite size 
effects. 

We analyzed the impact of the overpopulated high- 
energy tail on the cooling rate 7 which is the main hy- 
drodynamic coefficient of granular gases in the HCS. We 
found systematic deviations from the theoretical expres- 
sion which neglects the exponential tail. These devia- 
tions grow with increasing dissipation (decreasing e) as 
exp(— 6£), due to enhanced contributions from the tail. 

Finally, we observed and explained theoretically the 
extraordinary slow (as compared with molecular gases) 
relaxation of the velocity distribution to its asymptotic 
stationary form. It takes about ^ 20 — 30 collisions per 
particle and may be understood from the mechanism of 
the tail formation. 
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